using dnAnalytics.Math;
using NUnit.Framework;

namespace dnAnalytics.Tests.Math
{
    [TestFixture]
    public class SpecialFunctionsTests
    {
        private const double mAcceptableError = 1e-12;

        [Test]
        public void GammaLn()
        {
            Assert.AreEqual(System.Math.Log(1.0), SpecialFunctions.GammaLn(1.0), mAcceptableError);
            Assert.AreEqual(System.Math.Log(1.0), SpecialFunctions.GammaLn(2.0), mAcceptableError);
            Assert.AreEqual(System.Math.Log(2.0), SpecialFunctions.GammaLn(3.0), mAcceptableError);
            Assert.AreEqual(System.Math.Log(6.0), SpecialFunctions.GammaLn(4.0), mAcceptableError);
            Assert.AreEqual(System.Math.Log(24.0), SpecialFunctions.GammaLn(5.0), mAcceptableError);

            Assert.AreEqual(System.Math.Log(System.Math.Sqrt(System.Math.PI)), SpecialFunctions.GammaLn(0.5), mAcceptableError);
            Assert.AreEqual(System.Math.Log(System.Math.Sqrt(System.Math.PI) / 2.0), SpecialFunctions.GammaLn(1.5), mAcceptableError);
        }

        [Test]
        public void Gamma()
        {
            Assert.AreEqual(double.NaN, double.NaN);
            Assert.AreEqual(1.0, SpecialFunctions.Gamma(1.0), mAcceptableError);
            Assert.AreEqual(1.0, SpecialFunctions.Gamma(2.0), mAcceptableError);
            Assert.AreEqual(2.0, SpecialFunctions.Gamma(3.0), mAcceptableError);
            Assert.AreEqual(6.0, SpecialFunctions.Gamma(4.0), mAcceptableError);
            Assert.AreEqual(24.0, SpecialFunctions.Gamma(5.0), mAcceptableError);

            Assert.AreEqual(System.Math.Sqrt(System.Math.PI), SpecialFunctions.Gamma(0.5), mAcceptableError);
            Assert.AreEqual(System.Math.Sqrt(System.Math.PI) / 2.0, SpecialFunctions.Gamma(1.5), mAcceptableError);
        }

        [Test]
        public void DiGamma()
        {
            double EulerMascheroni = 0.57721566490153286;

            Assert.AreEqual(double.NaN, double.NaN);
            Assert.AreEqual(-EulerMascheroni, SpecialFunctions.DiGamma(1.0), mAcceptableError);
            Assert.AreEqual(-2.0 * System.Math.Log(2.0) - EulerMascheroni, SpecialFunctions.DiGamma(0.5), mAcceptableError);
            Assert.AreEqual(-System.Math.PI / 2.0 - 3.0 * System.Math.Log(2.0) - EulerMascheroni, SpecialFunctions.DiGamma(0.25), mAcceptableError);
        }

        [Test]
        public void BetaLn()
        {
            Assert.AreEqual(System.Math.Log(0.5), SpecialFunctions.BetaLn(1.0, 2.0), mAcceptableError);
            Assert.AreEqual(System.Math.Log(1.0), SpecialFunctions.BetaLn(1.0, 1.0), mAcceptableError);
        }

        [Test]
        public void Beta()
        {
            Assert.AreEqual(0.5, SpecialFunctions.Beta(1.0, 2.0), mAcceptableError);
            Assert.AreEqual(1.0, SpecialFunctions.Beta(1.0, 1.0), mAcceptableError);
        }

        [Test]
        public void Factorial()
        {
            Assert.AreEqual(1.0, SpecialFunctions.Factorial(1), mAcceptableError);
            Assert.AreEqual(2.0, SpecialFunctions.Factorial(2), mAcceptableError);
            Assert.AreEqual(6.0, SpecialFunctions.Factorial(3), mAcceptableError);
            Assert.AreEqual(24.0, SpecialFunctions.Factorial(4), mAcceptableError);
            Assert.AreEqual(120.0, SpecialFunctions.Factorial(5), mAcceptableError);
        }

        [Test]
        public void FactorialLn()
        {
            Assert.AreEqual(System.Math.Log(1.0), SpecialFunctions.FactorialLn(1), mAcceptableError);
            Assert.AreEqual(System.Math.Log(2.0), SpecialFunctions.FactorialLn(2), mAcceptableError);
            Assert.AreEqual(System.Math.Log(6.0), SpecialFunctions.FactorialLn(3), mAcceptableError);
            Assert.AreEqual(System.Math.Log(24.0), SpecialFunctions.FactorialLn(4), mAcceptableError);
            Assert.AreEqual(System.Math.Log(120.0), SpecialFunctions.FactorialLn(5), mAcceptableError);
        }

        [Test]
        public void Binomial()
        {
            Assert.AreEqual(1, SpecialFunctions.Binomial(1, 1), mAcceptableError);
            Assert.AreEqual(10, SpecialFunctions.Binomial(5, 2), mAcceptableError);
            Assert.AreEqual(35, SpecialFunctions.Binomial(7, 3), mAcceptableError);
        }

        [Test]
        public void Multinomial()
        {
            Assert.AreEqual(1, SpecialFunctions.Multinomial(1, new int[] { 1, 0 }), mAcceptableError);
            Assert.AreEqual(10, SpecialFunctions.Multinomial(5, new int[] { 3, 2 }), mAcceptableError);
            Assert.AreEqual(10, SpecialFunctions.Multinomial(5, new int[] { 2, 3 }), mAcceptableError);
            Assert.AreEqual(35, SpecialFunctions.Multinomial(7, new int[] { 3, 4 }), mAcceptableError);
        }      

        [Test]
        public void Erf()
        {
            Assert.AreEqual(1.0, SpecialFunctions.Erf(double.PositiveInfinity));
            Assert.AreEqual(-1, SpecialFunctions.Erf(double.NegativeInfinity));
            Assert.AreEqual(0.0, SpecialFunctions.Erf(0.0));
            TestHelper.TestSignificatDigits(1.0, SpecialFunctions.Erf(7.2), 15);
            TestHelper.TestSignificatDigits(-1.0, SpecialFunctions.Erf(-7.2), 15);
            TestHelper.TestSignificatDigits(0.842700792949715, SpecialFunctions.Erf(1.0), 15);
            TestHelper.TestSignificatDigits(0.520499877813047, SpecialFunctions.Erf(.5), 15);
            TestHelper.TestSignificatDigits(-0.842700792949715, SpecialFunctions.Erf(-1.0), 15);
            TestHelper.TestSignificatDigits(-0.520499877813047, SpecialFunctions.Erf(-.5), 15);

            TestHelper.TestSignificatDigits(1.12837916709551e-10, SpecialFunctions.Erf(1e-10), 15);
            TestHelper.TestSignificatDigits(0.112462916018285, SpecialFunctions.Erf(0.1), 15);
            TestHelper.TestSignificatDigits(0.222702589210478, SpecialFunctions.Erf(0.2), 15);
            TestHelper.TestSignificatDigits(0.328626759459127, SpecialFunctions.Erf(0.3), 15);
            TestHelper.TestSignificatDigits(0.428392355046668, SpecialFunctions.Erf(0.4), 15);
            TestHelper.TestSignificatDigits(0.520499877813047, SpecialFunctions.Erf(0.5), 15);
            TestHelper.TestSignificatDigits(0.603856090847926, SpecialFunctions.Erf(0.6), 15);
            TestHelper.TestSignificatDigits(0.677801193837418, SpecialFunctions.Erf(0.7), 15);
            TestHelper.TestSignificatDigits(0.742100964707661, SpecialFunctions.Erf(0.8), 15);
            TestHelper.TestSignificatDigits(0.796908212422832, SpecialFunctions.Erf(.9), 15);
            TestHelper.TestSignificatDigits(0.842700792949673, SpecialFunctions.Erf(0.9999999999999), 15);
            TestHelper.TestSignificatDigits(0.999991123853632, SpecialFunctions.Erf(System.Math.PI), 15);
            TestHelper.TestSignificatDigits(0.999999999999999, SpecialFunctions.Erf(5.7), 15);

            TestHelper.TestSignificatDigits(-1.12837916709551e-10, SpecialFunctions.Erf(-1e-10), 15);
            TestHelper.TestSignificatDigits(-0.112462916018285, SpecialFunctions.Erf(-0.1), 15);
            TestHelper.TestSignificatDigits(-0.222702589210478, SpecialFunctions.Erf(-0.2), 15);
            TestHelper.TestSignificatDigits(-0.328626759459127, SpecialFunctions.Erf(-0.3), 15);
            TestHelper.TestSignificatDigits(-0.428392355046668, SpecialFunctions.Erf(-0.4), 15);
            TestHelper.TestSignificatDigits(-0.520499877813047, SpecialFunctions.Erf(-0.5), 15);
            TestHelper.TestSignificatDigits(-0.603856090847926, SpecialFunctions.Erf(-0.6), 15);
            TestHelper.TestSignificatDigits(-0.677801193837418, SpecialFunctions.Erf(-0.7), 15);
            TestHelper.TestSignificatDigits(-0.742100964707661, SpecialFunctions.Erf(-0.8), 15);
            TestHelper.TestSignificatDigits(-0.796908212422832, SpecialFunctions.Erf(-.9), 15);
            TestHelper.TestSignificatDigits(-0.842696641833230, SpecialFunctions.Erf(-0.99999), 15);
            TestHelper.TestSignificatDigits(-0.999991123853632, SpecialFunctions.Erf(-System.Math.PI), 15);
            TestHelper.TestSignificatDigits(-0.999999999999999, SpecialFunctions.Erf(-5.7), 15);
        }

        [Test]
        public void Erfc()
        {
            Assert.AreEqual(0.0, SpecialFunctions.Erfc(double.PositiveInfinity));
            Assert.AreEqual(1.0, SpecialFunctions.Erfc(0));
            Assert.AreEqual(2.0, SpecialFunctions.Erfc(double.NegativeInfinity));
            Assert.AreEqual(2.0, SpecialFunctions.Erfc(-7.2), 1e-10);
            Assert.AreEqual(0.0, SpecialFunctions.Erfc(7.2), 1e-10);
            Assert.AreEqual(0.1572992071, SpecialFunctions.Erfc(1.0), 1e-10);
            Assert.AreEqual(0.4795001222, SpecialFunctions.Erfc(.5), 1e-10);
        }

        [Test]
        public void ErfInv()
        {
            Assert.AreEqual(double.NegativeInfinity, SpecialFunctions.ErfInv(-1.0));
            Assert.AreEqual(0.0, SpecialFunctions.ErfInv(0));
            Assert.AreEqual(double.PositiveInfinity, SpecialFunctions.ErfInv(1.0));
            TestHelper.TestSignificatDigits(8.86226925452758e-11, SpecialFunctions.ErfInv(1e-10), 15);
            TestHelper.TestSignificatDigits(0.0888559904942577, SpecialFunctions.ErfInv(0.1), 15);
            TestHelper.TestSignificatDigits(0.179143454621292, SpecialFunctions.ErfInv(0.2), 15);
            TestHelper.TestSignificatDigits(0.272462714726754, SpecialFunctions.ErfInv(0.3), 15);
            TestHelper.TestSignificatDigits(0.370807158593558, SpecialFunctions.ErfInv(0.4), 15);
            TestHelper.TestSignificatDigits(0.47693627620447, SpecialFunctions.ErfInv(0.5), 15);
            TestHelper.TestSignificatDigits(0.595116081449995, SpecialFunctions.ErfInv(0.6), 15);
            TestHelper.TestSignificatDigits(0.732869077959217, SpecialFunctions.ErfInv(0.7), 15);
            TestHelper.TestSignificatDigits(0.906193802436823, SpecialFunctions.ErfInv(0.8), 15);
            TestHelper.TestSignificatDigits(1.16308715367667, SpecialFunctions.ErfInv(.9), 15);
            TestHelper.TestSignificatDigits(3.4589107372755, SpecialFunctions.ErfInv(0.999999), 13);
            TestHelper.TestSignificatDigits(5.67591573974471, SpecialFunctions.ErfInv(.999999999999999), 13);


            TestHelper.TestSignificatDigits(-8.86226925452758e-11, SpecialFunctions.ErfInv(-1e-10), 15);
            TestHelper.TestSignificatDigits(-0.0888559904942577, SpecialFunctions.ErfInv(-0.1), 15);
            TestHelper.TestSignificatDigits(-0.179143454621292, SpecialFunctions.ErfInv(-0.2), 15);
            TestHelper.TestSignificatDigits(-0.272462714726754, SpecialFunctions.ErfInv(-0.3), 15);
            TestHelper.TestSignificatDigits(-0.370807158593558, SpecialFunctions.ErfInv(-0.4), 15);
            TestHelper.TestSignificatDigits(-0.47693627620447, SpecialFunctions.ErfInv(-0.5), 15);
            TestHelper.TestSignificatDigits(-0.595116081449995, SpecialFunctions.ErfInv(-0.6), 15);
            TestHelper.TestSignificatDigits(-0.732869077959217, SpecialFunctions.ErfInv(-0.7), 15);
            TestHelper.TestSignificatDigits(-0.906193802436823, SpecialFunctions.ErfInv(-0.8), 15);
            TestHelper.TestSignificatDigits(-1.16308715367667, SpecialFunctions.ErfInv(-.9), 15);
            TestHelper.TestSignificatDigits(-3.4589107372755, SpecialFunctions.ErfInv(-0.999999), 15);
            TestHelper.TestSignificatDigits(-5.67591573974471, SpecialFunctions.ErfInv(-.999999999999999), 15);
        }

        [Test]
        public void ErfcInv()
        {
            Assert.AreEqual(double.NegativeInfinity, SpecialFunctions.ErfcInv(2.0));
            Assert.AreEqual(double.PositiveInfinity, SpecialFunctions.ErfcInv(0.0));
            Assert.AreEqual(0.0, SpecialFunctions.ErfcInv(1));

            TestHelper.TestSignificatDigits(4.5728249673895, SpecialFunctions.ErfcInv(1e-10), 14);
            TestHelper.TestSignificatDigits(6.6015806223551, SpecialFunctions.ErfcInv(1e-20), 14);
            TestHelper.TestSignificatDigits(8.1486162231699, SpecialFunctions.ErfcInv(1e-30), 14);
            TestHelper.TestSignificatDigits(15.065574702593, SpecialFunctions.ErfcInv(1e-100), 13);
        }

        /// <summary>
        /// Compute the t'th harmonic number using a loop.
        /// </summary>
        private double ExactHarmonic(int t)
        {
            double r = 0.0;
            for (int i = 1; i <= t; i++)
            {
                r += 1.0 / i;
            }
            return r;
        }

        [Test]
        public void Harmonic()
        {
            TestHelper.TestSignificatDigits(ExactHarmonic(1), SpecialFunctions.Harmonic(1), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(2), SpecialFunctions.Harmonic(2), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(4), SpecialFunctions.Harmonic(4), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(8), SpecialFunctions.Harmonic(8), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(16), SpecialFunctions.Harmonic(16), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(100), SpecialFunctions.Harmonic(100), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(1000), SpecialFunctions.Harmonic(1000), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(10000), SpecialFunctions.Harmonic(10000), 14);
            TestHelper.TestSignificatDigits(ExactHarmonic(100000), SpecialFunctions.Harmonic(100000), 13);
            TestHelper.TestSignificatDigits(ExactHarmonic(1000000), SpecialFunctions.Harmonic(1000000), 12);
        }
    }
}